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We introduce a representation of electron operators as a product of a spin-carrying fermion and 
of a phase variable dual to the total charge (slave quantum rotor). Based on this representation, 
a new method is proposed for solving multi-orbital Anderson quantum impurity models at finite 
interaction strength U. It consists in a set of coupled integral equations for the auxiliary field Green's 
functions, which can be derived from a controlled saddle-point in the limit of a large number of field 
components. In contrast to some finite-C/ extensions of the non-crossing approximation, the new 
method provides a smooth interpolation between the atomic limit and the weak-coupling limit, and 
does not display violation of causality at low-frequency. We demonstrate that this impurity solver 
can be applied in the context of Dynamical Mean-Field Theory, at or close to half-filling. Good 
agreement with established results on the Mott transition is found, and large values of the orbital 
degeneracy can be investigated at low computational cost. 
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I. INTRODUCTION 



quantum impurity model (AIM) and its 
play a key role in several recent devel- 
field of strongly correlated electron sys- 
slectron devices, it has been widely used 
model for the competition between the 
de and the effect of tunnelling jij. In 
the Dynamical Mean-Field Theory 
gly correlated electron systems replaces 
system by an Anderson impurity 
;lf-consistently determined bath of con- 
Naturally, the AIM is also essen- 
rstanding of local moment formation in 
1 hat of heavy-fcrmion materials, particu- 
3d valence regime 0. 
important to have at our disposal quanti- 
g for the calculation of physical quan- 
with the AIM. The quantity of interest 
specific context. Many recent applica- 
^alculation of the localized level Green's 
:tral function), and possibly of some two- 
functions. 

impurity solvers" have been developed 
Broadly speaking, these methods fall in 
numerical algorithms which attempt at a 
on the computer, and analytical approx- 
i (which may also require some numeri- 
Among the former, the Quantum 
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Monte Carlo approach (based on the Hirsch-Fye algo- 
rithm) and the Numerical Renormalization Group play a 
prominent role. However, such methods become increas- 
ingly costly as the complexity of the model increases. In 
particular, the rapidly developing field of applications of 
DMFT to electronic structure calculations (j| requires 
methods that can handle impurity models involving or- 
bital degeneracy. For example, materials involving f- 
electrons are a real challenge to DMFT calculations when 
using Quantum Monte Carlo. Similarly, cluster exten- 
sions of DMFT require to handle a large number of cor- 
related local degrees of freedom. The dimension of the 
Hilbert space grows exponentially with the number of 
these local degrees of freedom, and this is a severe limi- 
tation to all methods, particularly exact diagonalization 
and numerical renormalization group. 

The development of fast and accurate (although ap- 
proximate) impurity solvers is therefore essential. In the 
one-orbital case.the iterated perturbation theory (IPT) 
approximation H has played a key role in the develop- 
ment of DMFT, particularly in elucidating the nature of 
the Mott transition |Q. A key reason for the success of 
this method is that it becomes exact both in the limit of 
weak interactions, and in the atomic limit. Unfortunately 
however, extensions of this approach to the multi-orbital 
case have not been as successful S. Another widely 
employed method is the non-crossing approximation, or 
NCA H pX[ . This method takes the atomic limit as a 
starting point, and performs a self-consistent resumma- 
tion of the perturbation theory in the hybridization to the 
conduction bath. It is thus intrinsically a strong-coupling 
approach, and indeed it is in the strong-coupling regime 
that the NCA has been most successfully applied. The 
NCA does suffer from some limitations however, which 
can become severe for some specific applications. These 
limitations are of two kinds: 

• The low-energy behaviour of NCA integral equa- 
tions is well-known to display non-Fermi liquid 
power laws. This can be better understood when 
formula ting the NCA approach in terms of slave 
bosons |]IT[ O] , It becomes clear then that the 
NCA actually describes accurately the overscreened 
regime of multichannel models. This is in a sense a 
remarkable success of NCA, but also calls for some 
care when applying the NCA to Fermi-liquid sys- 
tems in the screened regime. It must be noted that 
recent progresses have been made to improve the 
low-energy behaviour in the Fermi liquid case [fl3| . 

• A more important limitation for practical applica- 
tions has to do with the finite-?/ extensions of NCA 
equations. Standard extensions do not reproduce 
correctly the non- interacting U = limit. In con- 
trast to IPT, they are not "interpolative solvers" 
between the weak-coupling and the strong coupling 
limit. Furthermore, the physical self-energy tends 
to develop non-causal behaviour at low-frequency 
(i.e. X^'(^) > 0) below some temperature (the 



"NCA pathology") At half-filling and large U, this 
only happens at a rather low-energy scale, but away 
from half-filling or for smaller U, this scale can 
become comparable to the bandwith, making the 
finite- U NCA of limited applicability. 

In this article, we introduce new impurity solvers which 
overcome some of these difficulties. Our method is based 
on a rather general representation of strongly correlated 
electron systems, which has potential applications to lat- 
tice models as well . The general idea is to introduce a 
new slave-particle representation of physical electron op- 
erators, which emphasizes the phase variable dual to the 
total charge on the impurity. This should be contrasted to 
slave-boson approaches to a multi-orbital AIM: there, one 
introduces as many auxiliary bosons as there are fermion 
states in the local Hilbert space. This is far from eco- 
nomical: when the local interaction depends on the total 
charge only, it should be possible to identify a collective 
variable which provides a minimal set of collective slave 
fields. We propose that the phase dual to the total charge 
precisely plays this role. This turns a correlated electron 
model (at finite U) into a model of spin- (and orbital-) 
carrying fermions coupled to a quantum rotor degree of 
freedom. Various types of approximations can then be 
made on this model. In this article, we emphasize an ap- 
proximate treatment based on a sigma-model representa- 
tion of the rotor degree of freedom, which is then solved 
in the limit of a large number of components. This results 
in coupled integral equations which share some similar- 
ities to those of the NCA, but do provide the following 
improvements: 

• The non-interacting (U — 0) limit is reproduced ex- 
actly. For a fully symmetric multi-orbital model at 
half-filling and in the low temperature range, the 
atomic limit is also captured exactly, so that the 
proposed impurity solver is an interpolative scheme 
between weak and strong coupling. The sigma- 
model approximation does not treat as nicely the 
atomic limit far from half-filling howe ver (though 
improvements are possible, see section VE). 



• The physical self-energy does not display any viola- 
tion of causality (even at low temperature or small 
U). This is guaranteed by the fact that our inte- 
gral equations become exact in the formal limit of 
a large number of orbitals and components of the 
sigma-model field. 

It should be emphasized, though, that the low-energy 
behaviour of our equations is similar to the (infinite-?/) 
NCA, and characterized by non-Fermi liquid power laws 
below some low-energy scale. 

As a testing ground for the new solver, we apply it in 
this paper to the DMFT treatment of the Mott tran- 
sition in the multi-orbital Hubbard model. We find 
an overall very good agreement with the general as- 
pects of this problem, as known from numerical work 
and from some recently derived exact results. We also 
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compare to other impurity solvers (IPT, Exact Diago- 
nalization, QMC, NCA), particularly regarding the one- 
electron spectral function. 

This article is organized as follows. In section II, we 
introduce a representation of fermion operators in terms 
of the phase variable dual to the total charge (taking the 
finite U, multi-orbital Anderson impurity model as an 
example). In section III, a sigma- model representation 
of this phase variable is introduced, along with a gener- 
alization from 0(2) to 0(2M). In the limit of a large 
number of components, a set of coupled integral equa- 
tions is derived. In section |v|, we test this "dynamical 
slave rotor" (DSR) approach on the single-impurity An- 
derson model. The rest of the paper (Sec. |v|) is devoted 
to applications of this approach in the context of DMFT, 
which puts in perspective the advantages and limitations 
of the new method. 



methods jl4|, |l5), a bosonic field is introduced for each 
atomic state \<T\ ■ ■ ■ <tq) (along with spin-carrying auxil- 
iary fermions /J). Hence, these methods are not describ- 
ing the atomic spectrum in a very economical manner. 

The spectrum of (Q) can actually be reproduced by 
introducing, besides the set of auxiliary fermions fl, a 
single additional variable, namely the angular momen- 
tum L = —id /dO associated with a quantum 0(2) rotor 
9 {i.e. an angular variable in [0, 2-7r]. Indeed, the energy 
levels (§ can be obtained using the following hamiltonian 



local 



£0 Pa fa 



u 



L 



(3) 



A constraint must be imposed, which insures that the 
total number of fermions is equal to the 0(2) angular 
momentum (up to a shift): 



II. ROTORIZATION 

The present article emphasizes the role played by the 
total electron charge, and its conjugate phase variable. 
We introduce a representation of the physical electron 
in terms of two auxiliary fields: a fermion field which 
carries spin (and orbital) degrees of freedom, and the 
(total) charge raising and lowering operators which we 
represent in terms of a phase degree of freedom. The 
latter plays a role similar to a slave boson: here a "slave" 
0(2) quantum rotor is used rather than a conventional 
bosonic field. 



A. Atomic model 

We first explain this construction on the simple exam- 
ple of an atomic problem, consisting in an ./V-fold degen- 
erate atomic level subject to a local S'?7(iV)-symmetric 
Coulomb repulsion: 



H] 



local 



e d\d a + — 



(i) 



We use here eo = e<i + U/2 as a convenient redefinition 
of the impurity level, and we recast the spin and orbital 
degrees of freedom into a single index a — 1 ... N (for 
N = 2, we have a single orbital with two possible spin 
states a =],[)■ We note that eo is zero at half filling, 
due to particle-hole symmetry. 

The crucial point is that the spectrum of the atomic 
hamiltonian (|l|) depends only on the total fermionic 
charge Q = 0, • • • ,N and has a simple quadratic depen- 
dence on Q: 



Eq = e Q + 



U 



N 
~2 



(2) 



There are 2 N states, but only N + 1 different energy lev- 
els, with degeneracies (q) . In conventional slave boson 



L = 



E 



(4) 



This restricts the allowed values of the angular momen- 
tum to be I = Q - N/2 = -N/2, -N/2 + 1, • ■ • , N/2 - 
1, N/2, while in the absence of any constraint I can be an 
arbitrary (positive or negative) integer. The spectrum of 
(g) is e Q Q + Ul 2 /2, with I = Q - N/2 thanks to (§, so 
that it coincides with (g). 

To be complete, we must show that each state in the 
Hilbert space can be constructed in terms of these aux- 
iliary degrees of freedom, in a way compatible with the 
Pauli principle. This is achieved by the following identi- 
fication: 



\a 1 ...n Q ) d = \u x ...a Q ) f \S> = Q-N/2) 



(5) 



in which \<j\ . . . &q) d ^ denotes the antisymmetric fermion 
state built out of d— and /— fermions, respectively, and 
\£}g denotes the quantum rotor eigenstate with angular 
momentum /, i.e. (8 \i) g — e M . The creation of a physi- 
cal electron with spin a corresponds to acting on such a 
state with /J as well as raising the total charge (angular 
momentum) by one unit. Since the raising operator is 
e , this leads to the representation: 



d\ = lle M , d a = f a 



(6) 



Let us illustrate this for N — 2 by writing the four pos- 
sible states in the form: |t) d = \]) f |0) e , \[) d = ||) / |0) fl , 
\U) d = \U) f and |0) d = |0) / |-l) e , and showing 



that this structure is preserved by d\ — fie 10 . Indeed: 



IU)„ = 4li> I , = / T t |i) / e +tf |0) e = |U) / 



-1) 



(7) 



The key advantage of the quantum rotor representation 
is that the original quartic interaction between fermions 
has been replaced in (|^) by a simple kinetic term (L 2 ) 
for the phase field. 
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B. Rotor representation of Anderson impurity 
models 



C. Slave rotors, Hubbard-Stratonovich and gauge 
transformations 



We now turn to an SU (iV)-syrrrmetric Anderson im- 
purity model in which the atomic orbital is coupled to a 
conduction electron bath : 



U 



H = E e °4<k + j 



I 2 



(8) 



dl 



Cfc,, 



k.a 



Using the representation (g), we can rewrite this hamil- 
tonian in terms of the (/J, 9) fields only: 



U 



h = E e °/^ + ^ 2+ E e 4 



(9) 



k.a 



We then set up a functional integral formalism for the /J 
and degrees of freedom, and derive the action associated 
with (p~0|) - This is simply done by switching from phase 
and angular-momentum operators (9, L) to fields {9, d T 9) 
depending on imaginary time r 6 [0, /?]. The action is 
constructed from S = J^dr[—iL d T 9 + H + pd T f], and 
an integration over L is performed. It is also necessary 
to introduce a complex Lagrange multiplier h in order 
to implement the constraint L = J2 a — N/2. We 
note that, because of the charge conservation on the local 
impurity, h can be chosen to be independent of time, with 
ih e [0,271-/$. 

This leads to the following expression of the action: 



(d T 9 + ihf 



N, 



S = JdrY^fl{d T + e a -h)f a , >r 

+ E [ c lj d - + e ^ + V * c lj°e- i6 + h - c ] (10) 



We can recast this formula in a more compact form by 
introducing the hybridization function: 



A(iw) = Y - 



\V k \' 



iuj - e k 



fill 



and integrating out the conduction electron bath. This 
leads to the final form of the action of the SU (N) An- 
derson impurity model in terms of the auxiliary fermions 
and phase field: 



S = drJ2fl(dr + e a -h)f a + 





[dr[dr'A(T-T') V fl{T)Ur') e * (12) 
Jo Jo 



(d T 9 + ih) 2 N 

— — H h 

2U 2 



In this section, we present an alternative derivation of 
the expression ([l^) of the action which does not rely on 
the concept of slave particles. This has the merit to give a 
more explicit interpretation of the phase variable intro- 
duced above, by relating it to a Hubbard-Stratonovich 
decoupling field. This section is however not essential 
to the rest of the paper, and can be skipped upon first 
reading. 

Let us start with the imaginary time action of the An- 
derson impurity model in terms of the physical electron 
field for the impurity orbital: 



S = \ rfrV;4(9 T + eoK 
Jo 



U 
~2 



E4*t" 



N 



[dr[dr'A(r-r')J24(r)dAr') 
Jo Jo 



(13) 



Because we have chosen a SU (N )-symmetric form for the 
Coulomb interaction, we can decouple it with only one 
bosonic Hubbard-Stratonovitch field 6(t): 



S = [drJ24(dr 
Jo 



+ e + i4>{ T ))d a + 



fir) .N 
2U 



[dr[dr'A(T-T')J24(T)d a (T') 
Jo Jo 



(14) 



Hence, a linear coupling of the field 4>(t) to the fermions 
has been introduced. The idea is now to eliminate 
this linear coupling for all the Fourier modes of the </>- 
field, except that corresponding to zero-frequency: (f>o = 
Jq 4> [27r] . This can be achieved by performing the fol- 
lowing gauge transformation: 



4( T ) = /t( T ) e *Jo r ^ e - l, W' 3 



(15) 



The reason for the second phase factor in this expression 
is that it guarantees that the new fermion field /J also 
obeys antiperiodic boundary conditions in the path in- 
tegral. It is easy to check that this change of variables 
does not provide any Jacobian, so that the action simply 
reads: 



10 

rP r P 



o w , 4> 2 {r) N<f> Q 
2U 1 2 [3 



[ drl dr'A(r-r') V/](r)/ CT (r')e l ^^-^l (16) 
Jo Jo 



We now set: 



d9 1 f p \ 

M = + p ^ I with: ^ = J [2tt] J (17) 
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and notice that the field 6(t) has the boundary condition 
0(0) = 0(0) [2tt]. It therefore corresponds to an 0(2) 
quantum rotor, and the expression of the action finally 
reads: 







S = /dr^/t(9 T + e 



o 

■ [dr [dr'A( T -T') T fl(r)Ur')e^-^' 
Jo Jo 



• tow , (^fl+f) 2 N^q 

1 I3 )la+ 2U 1 2 H 



(18) 



This is exactly expression (|I2[), with the identification: 
4>o/(3 = i/i- This, together with fll^), provides an explicit 
relation between the quantum rotor and Lagrange mul- 
tiplier fields on one side, and the Hubbard-Stratonovich 
field conjugate to the total charge, on the other. 



action reads 



S= drY i ti{a r + e -h)f <r + —h-M 
Jo 

Jo 



2U 



MX 



\d T X a \ 
2U 



— (X*d T X a -h.c.) + X\X a \ 2 



+ fdr[dT'^A(r-T')J2ftir)fAr')X a (T)X* a (r') 



Jo 



Let us note that this action corresponds to the hamilto- 
nian: 



H 



2M 



a,/3 



k,aa 



2 ( C k,<rJ<rK + ft C k,<T«X a 



(21) 



III. SIGMA-MODEL REPRESENTATION AND 
SOLUTION IN THE LIMIT OF MANY 
COMPONENTS 



A. From quantum rotors to a sigma model 

Instead of using a phase field to represent the 0(2) 
degree of freedom, one can use a constrained (complex) 



bosonic field X 



r.iO 



with: 



\X(r) 



(19) 



The action ( |12[ ) can be rewritten in terms of this field, 
provided a Lagrange multiplier field A(r) is used to im- 
plement this constraint: 



fP AT U2 

S = J dT E + £ o - h )f° + Y h ~2U 

[dT[dr'A(T-r')J2fUr)fa(r')X(T)X*(r') (20) 
Jo Jo 



Hence, the Anderson model has been written as a theory 
of auxiliary fermions coupled to a non-linear 0(2) sigma- 
model, with a constraint (implemented by h) relating the 
fermions and the sigma- model field X(t). 

A widely-used limit in which sigma models become 
solvable, is the limit of a large number of components 
of the field. This motivates us to generalize (20) to a 
model with an 0(2M) symmetry. The bosonic field X 
is thus extended to an M-component complex field X a 
(a = 1 . . . M) with J2 a \X a \ 2 — M. The corresponding 



In this expression, L Q ,/3 denotes the angular momen- 
tum tensor associated with the X a vector. The hamilto- 
nian (fl]) is a generalization of the SU(N) x 0(2) = U(N) 
Anderson impurity model to an SU(N) x 0(2M) model 
in which the total electronic charge is associated with a 
specific component of L. It reduces to the usual Ander- 
son model for M = 1. 

In the following, we consider the limit where both N 
and M become large, while keeping a fixed ratio N/M . 
We shall demonstrate that exact coupled integral equa- 
tions can be derived in this limit, which determine the 
Green's functions of the fcrmionic and sigma-model fields 
(and the physical electron Green's function as well). The 
fact that these coupled integral equations do correspond 
to the exact solution of a well-defined hamiltonian model 
(Eq. (pl|)) guarantees that no unphysical features (like e.g 
violation of causality) arise in the solution. Naturally, the 
generalized hamiltonian (|2l]) is a formal extension of the 
Anderson impurity model of physical interest. Extending 
the charge symmetry from 0(2) to 0(2M) is not entirely 
inocuous, even at the atomic level: as we shall see below, 
the energy levels of a single <3(2M) quantum rotor have 
multiple degeneracies, and depend on the charge (angular 
momentum) quantum number in a way which does not 
faithfully mimic the 0(2) case. Nevertheless, the basic 
features defining the generalized model (a localized or- 
bital subject to a Coulomb charging energy and coupled 
to an electron bath by hybridization) are similar to the 
original model of physical interest. 



B. Integral equations 

In this section, we derive coupled integral equations 
which become exact in the limit where both M and N 
are large with a fixed ratio: 



(N, M -► oo) 



(22) 
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Following |12| (see also |Tl|), the two body interaction be- 
tween the auxiliary fermions and the sigma-modcl fields 
is decoupled using (bosonic) bi-local fields Q(t, t') and 
Q(t, t') depending on two times. Hence, we consider the 
action: 



P N u2 

dr V fUdr + eo - h)f a +—h- M^ttt - MX 



2U 







+ fa fa M Q(Ty)Q(r r> } 
Jo Jo A ( T - T ) 

rP rP 



fdr fdr' Q(r,r')^X a (r)X*(r') 
Jo Jo a 

fdr fdr'Q(r,r')J2fUr)fa(r') 
Jo Jo 



(23) 



In the limit fl22j), this action is controlled by a saddle- 
point, at which the Lagrange multipliers take static ex- 
pectation values h and A, while the saddle point values 
of the Q and Q fields are translation invariant functions 
of time Q(t — t') and Q(t — r'). 

Introducing the imaginary-time Green's functions of 
the auxiliary fermion and sigma-model fields as: 



G f (r) ee -(T T / CT (r)/t(0)) 
G x (r) = + (T T X a (r)X*(0)) 



(24) 
(25) 



we differentiate the effective action ( |23| ) with respect 
to Q(t) and Q(t), which leads to the following saddle- 
point equations: Q{t) — — 7V"A(r)G/(— r) and Q(t) = 
A(t)Gx(t). The functions Q(iu) n ) (= Ey) and Q (iv n ) 
(= Ex) define fermionic and bosonic self-energies: 

G^iicjn) = iu} n — e + h-T,f(iu! n ) (26) 

G~x{iVn) = £ + \-^-Vx(iu n ) (27) 

where oj n (resp. v n ) is a fermionic (resp. bosonic) Mat- 
subara frequency. The saddle point equations read: 



Sx(r) = -AfA(-r)G f (r) 
E ; (r) = A(r)G x (r) 



(28) 
(29) 



together with the constraints associated with h and A: 
G x (t = 0) = 1 (30) 

(31) 



+ j^j[d T G x (T = 0-) +d T G x (T = 0+)] 

There is a clear similarity between the structure of these 
coupled integral equations and the infinite-C/ NCA equa- 
tions J|. We note also significant differences, such as the 
constraint equations. Furthermore, the finite value of the 



Coulomb repulsion U enters the bosonic propagator (]27 
in a quite novel manner. 

The two key ingredients on which the present method 
are based is the use of a slave rotor representation of 
fermion operators, and the use of integral equations for 
the frequency-dependent self-energies and Green's func- 
tions. For this reason, we shall denote the integral equa- 
tions above under the name of "Dynamical Slave Rotor" 
method (DSR) in the following. 



C. Some remarks 

We make here some technical remarks concerning these 
integral equations. 

First, we clarify how the interaction parameter U was 
scaled in order to obtain the DSR equations above. This 
issue is related to the manner in which the atomic limit 
(A = 0) is treated in this method. In the original 
0(2) atomic hamiltonian ([!]), the charge gap between 
the ground-state and the first excited state is U/2 at 
half- filling (eo = 0). In the DSR method, the charge gap 
is associated with the gap in the slave rotor spectrum. If 
the 0(2M) generalization of (Q) is written as in (pl|): 



U 

2M 



a,/3 



(32) 



the spectrum reads: E e = U£(£ + 2M - 2)/(2M). As a 
result, the energy difference from the ground state to the 
first excited state is E 1 - E = U (2M - 1)/(2M) ~ U 
at large M, whereas it is U/2 at M = 1. In order to use 
the DSR method in practice as an approximate impurity 
solver, the parameter U should thus be normalized in a 
different way than in (H|), so that the gap is kept equal 
to U/2 in the large-M limit as well. Technically this can 
be enforced by choosing the following normalization: 



TTl _ 

n int — 



u 



AM 



1 

a,0 



(33) 



instead of (|32j ) . Note that this scaling coincides with ([32 
for M = 1, but does yield E 1 - E = U/2 for large-M, 
as desired. This definition of U was actually used when 
writing the saddle-point integral equations ( p7[ ) , although 
we postponed the discussion of this point to the present 
section for reasons of simplicity. 

Let us elaborate further on the accuracy of the DSR 
integral equations in the atomic limit. In Appendix 
we show that the physical electron spectral function ob- 
tained within DSR in the atomic limit coincides with the 
exact 0(2) result at half-filling and at T — 0. This is 
a non-trivial result, given the fact that the constraint is 
treated on average and the above remark on the spectrum 
spectrum. In contrast to NCA, the DSR method (in its 
present form) is not based by construction on a strong- 
coupling expansion around the exact atomic spectrum, 
so that this is a crucial check for the applicability of this 
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method in practice. In the context of DMFT for exam- 
ple, it is essential in order to describe correctly the Mott 
insulating state 0. However, the DSR integral equations 
fail to reproduce exactly the 0(2)atomic limit off half- 
filling, as explained in Appendix EL Deviations be come 
severe for too high dopings, as discussed in Sec. VE . 
This makes the present form of DSR applicable only for 
systems in the vicinity of half-filling. 

We now discuss some general spectral properties of the 
DSR solver. From the representation of the physical elec- 
tron field dl = fXX, and from the convention chosen for 
the pseudo-particules Green's functions (j24|-p5|), the one- 
electron physical Green's function is simply expressed as: 



G d (T)=G f (T)Gx(-r) 



(34) 



Therefore eq. (|30|), combined with the fact that /t has 
a (— 1) discontinuity at r = (which is obvious from 
(Pq)), shows that a% possesses also a (—1) jump at zero 
imaginary time. This ensures that the physical spec- 
tral weight is unity in our theory, and thus that phys- 
ical spectral function are correctly normalized. Because 
the DSR integral equations result from a controlled large 
N,M limit, it also insures that the physical self-energy 
always have the correct sign (i.e. ImY,d{u + i0 + ) < 0). 
This is not the case |H| for the finite U version of the 
NCA [|l0] which is constructed as a resummation of the 
strong-coupling expansion in the hybridization function 
A(r) (Strong coupling resummations generically suffer 
from non-causality, see also |l7j for an illustration.) 

Finally, we comment on the non-interacting limit U — ► 
0. This is a major failure of the usual NCA, which lim- 
its its applicability in the weakly correlated regime. In 
the DSR formalism, this limit is exact as can be noticed 
from equation (p7j). Indeed, as U vanishes, only the zero- 
frequency component of Gx(iv n ) survives, so that Gx{t) 
simply becomes a constant. Because of the constraint 
©, we get correctly G x (t) = 1 at U = 0. From @ 
and (HH), this proves that Gd(r) is the non- interacting 
Green's function: 



1 



- eo - A(iu„ 



(35) 



We finally acknowledge that an alternative dynamical 
approximation to the finite U Anderson model jl3| was 
recently developed as an extension of NCA by Kroha, 
Wolfle and collaborators (a conventional slave boson rep- 
resentation was used in this work). Many progresses have 
been made following this method, but, to the authors' 
knowledge, this technique has not yet been implemented 
in the context of DMFT (one of the reasons is its com- 
putational cost). By developing the DSR approximation, 
we pursue a rather complementary goal: the aim here is 
not to improve the low-energy singularities usually en- 
countered with integral equations, but rather to have a 
fast and efficient solver which reproduces correctly the 
main features of the spectral functions and interpolates 
between weak and strong coupling. In that sense, it is 
very well adapted to the DMFT context. 



IV. APPLICATION TO THE 
SINGLE-IMPURITY ANDERSON MODEL 

We now discuss the application of the DSR in the sim- 
plest setting: that of a single impurity hybridized with a 
fixed bath of conduction electrons. For simplicity, we fo- 
cus on the half-filled, particle-hole symmetric case, which 
implies £q = h = 0. The doped (or mixed valence) case 
will be addressed in the next section, in the context of 
DMFT. 

As the strength of the Coulomb interaction U is in- 
creased from weak to strong coupling, two well-known 
effects are expected (see e.g j|). First the width of the 
low-energy resonance is reduced from its non-interacting 
value Ao = |A"(0)|. As one enters the Kondo regime 
A <C U <C A (with A the conduction electrons band- 
width), this width becomes a very small energy scale, of 
the order of the Kondo temperature: 



T K = V2[/AoCxp(-7rl7(8Ao)) 



(36) 



A (local) Fermi liquid description applies, with quasi- 
particles having a large effective mass and small weight: 
Z = m/m* ~ Tk/ 'Ao. The impurity spin is screened for 
T < T K . 

Second, the corresponding spectral weight is trans- 
ferred to high energies, into "Hubbard bands" associated 
to the atomic-like transitions (adding or removing an 
electron into the half- filled impurity orbital), broadened 
by the hybridization to the conduction electron bath. 
The suppression of the low-energy spectral weight cor- 
responds to the suppression of the charge fluctuations on 
the local orbital. These satellites are already visible at 
moderate values of the coupling U/Aq. As temperature 
is increased from T < Tk to T > Tk, the Kondo quasi- 
particle resonance is quickly destroyed, and the missing 
spectral weight is added to the Hubbard bands. 

The aim of this section is to investigate whether the in- 
tegral equations introduced in this paper reproduce these 
physical effects in a satisfactory manner. 



A. Spectral functions 

We have solved numerically these integral equations by 
iteration, both on the imaginary axis and for real frequen- 
cies. Working on the imaginary axis is technically much 
easier. A discretization of the interval r € [0,/3] is used 
(with typically 8192 points, and up to 32768 for reaching 
the lowest temperatures), as well as Fast Fourier Trans- 
forms for the Green's functions. Searching by dichotomy 
for the saddle-point value of the Lagrange multiplier A 
(Eq. ^) is conveniently implemented at each step of the 
iterative procedure. Technical details about the analytic 
continuation of ( p8| - pl] ) to real frequencies, as well as their 
numerical solution are given in Appendix |^. 

On Fig. [I], we display our results for the impurity- 
orbital spectral function pd(oj), at three different temper- 
atures (the density of states of the conduction electron 
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FIG. 1: d-level spectral function Pd(w) for fixed U — 2 and in- 
verse temperatures j3 = 4, 20, 2000, with the conduction elec- 
tron bath described in the text 



0.06 




FIG. 2: Kondo temperature Tk from the exact formula (line) 
and from the numerical solution of the DSR equations (dots) . 
Units of energy are such that Ao = 0.16. 



bath is chosen as a semi-circle with half- width A = 6, the 
resonant level width is Ao = 0.16 and we take U = 2). 
The growth of the Kondo resonance as the temperature is 
lowered is clearly seen. The temperatures in Fig. |l| have 
been chosen such as to illustrate three different regimes: 
for T ^> Tr , no resonance is seen and the spectral density 
displays a "pseudogap" separating the two high-energy 
bands; for T ~ Tk transfer of spectral weight to low 
energy is seen, resulting in a fully developed Kondo res- 
onance for T <C Tk- We have not obtained an analyti- 
cal determination of the Kondo temperature within the 
present scheme. In the case of NCA equations, it is pos- 
sible to derive a set of differential equations in the limit 
of infinite bandwith which greatly facilitate this. This 
procedure cannot be applied here, because of the form 
( p7j ) of the boson propagator. Nevertheless, we checked 
that the numerical estimates of the width of the Kondo 
peak is indeed exponentially small in U as in formula 
(36), see Fig. |^. However, because U is normalized as in 
(33), (which gives the correct atomic limit), the prefactor 
inside the exponential appears to be twice too small. 

On Fig. ||, we display the spectral function for a fixed 
low temperature and increasing values of U . The strong 
reduction of the Kondo scale (resonance width) upon in- 
creasing U is clear on this figure. We note that the high- 
energy peaks have a width which remains of order Ao, 
independently of U, which is satisfactory. However, we 
also note that they are not peaked exactly at the atomic 
value ±Z7 /2, which might be an artefact of these integral 
equations. The shift is rather small however. 

Fig. H illustrates how the high- and low- energy fea- 
tures in the d-level spectral functions are associated with 
corresponding features in the auxiliary particle spectral 
functions pf and px- In particular, the sigma- model 
boson (slave rotor) is entirely responsible for the Hub- 
bard bands at high energy (as expected, since it describes 
charge fluctuations). 

As stressed in the introduction, an advantage of our 




FIG. 3: Pd(w) at low temperature (j3 = 600) and for increasing 
[7=1,2,3 




FIG. 4: Pseudo-particule spectral function at U = 2 and 
jfl = 600. The Kondo resonance is visible in pf(ui), broken 
curve, whereas px(u>) displays higher energy features, plain 
curve 
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scheme is that the d-level self-energy is always causal, 
even for small U or large doping. This is definitely an 
improvement as compared to the usual {/-NCA approx- 
imation. This is illustrated by Fig. [5], from which it is 
also clear that decreases (and eventually vanishes) as 
U goes to zero (for a more detailed discussion and com- 
parison to NCA, see Sec. HI). However, it is also clear 
from this figure that the low-energy behavior of the self- 
energy is not consistent with Fermi liquid theory. This 
is a generic drawback of NCA-like integral equation ap- 
proaches, that we now discuss in more details. 




FIG. 5: Imaginary part of the physical self-energy on the real- 
frequency axis, for U = 0.5 (upper curve) and U = 2 (lower 
curve) at @ = 600 



B. Low-energy behaviour and Friedel sum-rule 

We discuss here the low-energy behavior of the integral 
equations (in the case of a featureless conduction electron 
bath), for both the d-level and auxiliary field Green's 
functions. As explained below, this low-energy behavior 
depends sensitively on the ratio Af = N/M which is kept 
fixed in the limit considered in this paper. The calcula- 
tions are detailed in Appendix where we establish the 
following. 

• Friedel sum rule. The zero-frequency value of the 
d-level spectral function at T — is independent of 
U and reads: 



-1 7i72 /vr Af 

p d (uj = T = 0) = . / : tan [ -- 



ttA"(0)AA+1 \2Af+l 



(37) 



This is to be contrasted with the exact value for 
the 0(2) model (M = 1), which is independent of 
N and reads: 

/9r ct (w = r=0) = _l_ (M=1; anyAr) (38) 



This is the Friedel sum rule j|, p|, which in this 
particle-hole symmetric case simply follows from 
the Fermi-liquid requirement that the (inverse) life- 
time S^'(o; = 0) should vanish at T = (since 
G^ 1 = iuj — A(iuj) — Y*d(iuj)). As a result, pd(0) 
is pinned at its non-interacting value. The inte- 
gral equations discussed here yield a non- vanishing 
E^'(cj = 0) (albeit always negative in order to sat- 
isfy causality), and hence do not describe a Fermi 
liquid at low energy. The result ( |37j ) is identical to 
that found in the NCA for U = oo, but holds here 
for arbitrary U. There is actually no contradiction 
between this remark and the fact that our integral 
equations yield the exact spectral function in the 
t7 — ^ limit. Indeed, the limit U — ► at finite 
T, lo does not commute with oj, T — > at finite U 
(in which © holds). 

• Low-frequency behavior. The auxiliary particle 
spectral functions have a low-frequency singu- 
larity characterized by exponents which depend 
continuously on Af (as in the U — oo NCA): 
p f (w) cx l/E'/M ex l/\u\ a f, p x {w) cx 1/E£.(w) cx 
Sign(w)/M Q *, with: a f = 1 - a x = l/(Af + !)• 
These behavior are characterized more precisely in 
Appendix |c[ A power-law behavior is also found 
for the physical self-energy S^(tj) — S^'(0) at low 
frequency (as evident from Fig. ^). 

Let us comment on the origin of these low-energy fea- 
tures, as well as on their consequences for the practical 
use of the present method. 

First, it is clear from expression ( |2l| ) that the Ander- 
son impurity hamiltonian generalized to SU(7V) x 0(2M ) 
actually involves M channels of conduction electrons. 
Hence, the non-Fermi liquid behaviour found when solv- 
ing the integral equations associated with the N, M — > oo 
limit simply follows from the fact that multi-channel 
models lead to overscreening of the impurity spin, and 
correspond to a non- Fermi liquid fixed point. In that 
sense, these integral equations reproduce very accurately 
the expected low-energy physics, as previously studied 
for the simplest case of the Kondo model in jll], [l^] . 

Naturally, this means that the use of such integral 
equations to describe the one-channel (exactly screened) 
case becomes problematic in the low-energy region. In 
particular, the exact Friedel sum rule is violated, the d- 
level lifetime remains finite at low energy and non-Fermi 
liquid singularities are found. While the approach is rea- 
sonable in order to reproduce the overall features of the 
one-electron spectra, it should not be employed to calcu- 
late transport properties at low energy for example. We 
note however that the deviation from the exact Friedel 
sum rule vanishes in the Af — > oo limit. This is expected 
from the fact that in this limit the number of channels 
(M) is small as compared to orbital degeneracy (N). The 
violation of the sum rule remains rather small even for 
reasonable values of Af. This parameter can actually be 
used as an adjustable parameter when using the present 
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method as an approximate impurity solver. There is no 
fundamental reason for which JV = N should provide 
the best approximate description of the spectral func- 
tions of the one-channel case. We shall use this possi- 
bility when applying this method in the DMFT context 
in the next section: there, we choose M = 3 in order to 
adjust to the known critical value of U for a single or- 
bital. We also used this value in the calculations reported 
above. The zero-frequency value of the spectral density 
is thus Pd(0) — 1.85, while the Friedel sum rule would 
yield p<i(0) ~ 1.95 (cf. Fig. ^). Hence the violation of the 
sum rule is a small effect (of the order of 5%), comparable 
to the one found with "numerically exact" solvers, due 
to discretization errors |l9|] . Also, we point out that the 
pinning of pd(0) at a value independent of U (albeit not 
that of the Friedel sum-rule) is an important aspect of 
the present method, which will prove to be crucial in the 
context of DMFT in order to recover the correct scenario 
for the Mott transition. 

Finally, we emphasize that increasing the parameter 
M also corresponds to increasing the orbital degeneracy 
of the i mpur ity level. This will be studied in more details 
in Sec. |VD| . In particular, we shall see that correlation 
effects become weaker as M is increased (for a given value 
of U), due to enhanced orbital fluctuations. 



V. APPLICATIONS TO DYNAMICAL 
MEAN-FIELD THEORY AND THE MOTT 
TRANSITION 

A. One-orbital case: Mott transition, phase 
diagram 

Dynamical Mean-Field Theory has led to significant 
progress in our understanding of the physics of a corre- 
lated metal close to the Mott transition pj . The detailed 
description of this transition itself within DMFT is now 
well established @, @, |l|, |l|, |(], ||, ||. In this section, 
we use these established results as a benchmark and test 
the applicability of the method introduced in this paper 
in the context of DMFT, with very encouraging results. 
As explained above, this is particularly relevant in view 
of the recent applications of DMFT to electronic struc- 
ture calculations of correlated solids ||, which call for 
efficient multi-orbital impurity solvers. 

As is well-known 0, |L 0] , DMFT maps a lattice hamil- 
tonian onto a self-consistent quantum impurity model. 
We discuss first the half-filled Hubbard model, and ad- 
dress later the doped case. We then have to solve a 
particle-hole symmetric Anderson impurity model : 



subject to the self-consistency condition: 
A(r) = t 2 G d (r) 



(40) 



f/3 

Jo 



E 4<k - 1 



(39) 



In this expression, a semi-circular density of states with 
half-bandwith D = 2t has been considered, correspond- 
ing to a infinite-connectivity Bethe lattice (z — oo) with 
hopping tij — tj\fz. In the following, we shall generally 
express all energies in units of D (D — 1). 

In practice, one must iterate numerically the "DMFT 
loop": A(r) -> G d (r) -> A(r) new = t 2 G d , using some 
"impurity solver" . Here, we make use of the integral 
equations ( p8f]3l| ). The hybridization function A(r) be- 
ing determined by the self-consistency condition (licj), 
there are only two free parameters, the local Coulomb 
repulsion U and the temperature T (normalized by D). 

We display in Fig. || the spectral functions obtained at 
low temperature, for increasing values of U, and in Fig. |?j 
the corresponding phase diagram. The value of the pa- 
rameter Af has been adapted to the description of the 
one-orbital case (see below). The most important point 
is that we find a coexistence region at low-enough tem- 
perature: for a range of couplings U C \{T) <U< U C 2(T), 
both a metallic solution and an insulating solution of the 
(paramagnetic) DMFT equations exist. The Mott tran- 
sition is thus first-order at finite temperatures. This is in 
agreement with the results established for this problem 
by solving the DMFT equations with controlled numer- 
ical methods jl9[ as well as with analytical results 
p0| , p2[ . In particular, the spectral functions that we ob- 
tain (Fig. H) display the well-known separation of energy 
scales found within DMFT: there is a gradual narrowing 
of the quasiparticle peak, together with a preformed Mott 
gap at the transition. In the next section, we compare 
these spectral functions to those obtained using other 
approximate solvers. As pointed out there, despite some 
formal similarity in the method, it is well-known that 
the standard U-NCA does not reproduce correctly this 
separation of energy scales close to the transition [E3[ . 




fdr fdr' A(r-r')J24(r)d a (r') 
Jo Jo 



FIG. 6: Local spectral function at (3 = 40 and U = 1,2,3 
for the half-filled Hubbard model within DMFT, as obtained 
with the DSR solver 



11 



0.1 



0.075 



T 0.05 



0.025 



1 1 1 


1 
























\j \ 










1 , 1 





2 

u 



FIG. 7: Single-orbital (paramagnetic) phase diagram at half- 
filling. Squares indicate the ED result, the DSR result is the 
solid line, and IPT the broken line 



present method to those obtained with IPT and ED is 
displayed in Fig. ^, for a value of U corresponding to a 
correlated metal close to the transition. 




Let us explain how the parameter Af has been chosen 
in these calculations. As demonstrated below, the val- 
ues of the critical couplings U c i and U C 2 (and hence the 
whole phase diagram and coexistence window) strongly 
depends on the value of this parameter. This is ex- 
pected, since TV" is a measure of orbital degeneracy. What 
has been done in the calculations displayed above is 
to choose Af in such a way that the known value pf| 
U C 2{T — 0) ~ 2.9 of the critical coupling at which the 
T = metallic solution disappears in the single-orbital 
case, is accurately reproduced. We found that this re- 
quires Af ~ 3 (note that N/M = 2 in the one-orbital 
case, so that the best agreement is not found by a naive 
application of the large N,M limit). This value being 
fixed, we find a critical coupling U C \{T = 0) ~ 2.3 in 
good agreement with the value from (adaptative) exact 
diagonalizations U c \ ~ 2.4. The whole domain of coexis- 
tence in the (U, T) plane is also in good agreement with 
established results (in particular we find the critical end- 
point at T c ~ 1/30, while QMC yields T c ~ 1/40). These 
are very stringent tests of the applicability of the present 
method, since we have allowed ourselves to use only one 
adjustable parameter (Af). In Sec. VD, we study how the 
Mott transition depends on the number of orbitals, which 
further validates the procedure followed here. 



B. Comparison to other impurity solvers: spectral 
functions 



Let us now compare the spectral functions obtained 
by the present method with other impurity solvers com- 
monly used for solving the DMFT equations. We start 
with the iterated perturbation theory approximation 
(IPT) and the exact diagonalization method (ED). Both 
methods have played a major role in the early develop- 
ments of the DMFT approach to the Mott transition . 
A comparison of the spectral functions obtained by the 



FIG. 8: Comparison between DSR (dash), IPT (straight) and 
ED (dot) at U = 2.4 and (3 = 60 

The overall shape and characteristic features of the 
spectral function are quite similar for the three meth- 
ods. A narrow quasiparticle peak is formed, together 
with Hubbard bands, and there is a clear separation of 
energy scales between the width of the central peak (re- 
lated to the quasiparticle weight) and the "preformed" 
Mott (pseudo-)gap associated with the Hubbard bands: 
this is a distinguishing aspect of DMFT. It is crucial for 
an approximate solver to reproduce this separation of en- 
ergy scales in order to yield a correct description of the 
Mott transition and phase diagram. 

There are of course some differences between the three 
methods, on which we now comment. First, we note that 
the IPT approximation has a somewhat larger quasipar- 
ticle bandwith. This is because the transition point U C 2 
is overestimated within IPT (cf. Fig. 0), so that a more 
fair comparison should perhaps be made at fixed U/U C 2- 
It is true however, that the DSR method has a tendency 
to underestimate the quasiparticle bandwith, and par- 
ticularly at smaller values of U. Accordingly, the Hub- 
bard bands have a somewhat too large spectral weight, 
but are correctly located in first approximation. The 
detailed shape of the Hubbard bands is not very accu- 
rately known, in any case. (The ED method involves 
a broadening of the delta- function peaks obtained by 
diagonalizing the impurity hamiltonian with a limited 
number of effective orbitals, so that the high-energy be- 
haviour is not very accurate on the real axis. This is also 
true, actually, of the more sophisticated numerical renor- 
malization group). We emphasize that, since the DSR 
method does not have the correct low-frequency Fermi 
liquid behaviour, the quasiparticle bandwith should be 
interpreted as the width of the central peak in pd(w) 
(while the quasiparticle weight Z cannot be defined for- 
mally) . 

In Fig U and Fig. |l^, we display the spectral functions 
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FIG. 9: (7-NCA spectral function for U = 1.5 at low temper- 
ature. 




FIG. 10: Results for U = 0.1 showing how (7-NCA overshoots 
the Friedel's sum rule. The Slave Rotor method is correctly 
converging towards the free density of states (semi-circular) 



obtained by using the NCA method (extended at finite 
U in the simplest manner). It is clear that the [/-NCA 
underestimates considerably the quasiparticle bandwith 
(and thus yields a Mott transition at a rather low value 
of the coupling). This is not very surprising, since this 
method is based on a strong coupling expansion around 
the atomic limit, and one could think again of a com- 
parison for fixed U/U c . More importantly, the {7-NCA 
misses the important separation of energy scales between 
the central peak and the Mott gap. As a result, it does 
not reproduce correctly the phase diagram for the Mott 
transition within DMFT (in particular regarding coex- 
istence). A related key observation is that the [/-NCA 
does not have the correct weak-coupling (small U) limit. 
To illustrate this point, we display in Fig. |l^ the {/-NCA 
and DSR spectral functions for a tiny value of the inter- 
action U/D — 0.1: the DSR result is shown to approach 
correctly the semi-circular shape of the non-interacting 
density of states, while the [/-NCA displays a charac- 
teristic inverted V-shape: in this regime the violation 
of the Friedel sum rule becomes large and the negative 
lifetime pathology is encountered within {/-NCA. In con- 
trast, the DSR yields a pinning of the spectral density 
at a [/-independent value and does not lead to a viola- 
tion of causality. It should be emphasized however that, 
even though it yields the exact {7 = density of states 
at T — 0, the DSR method is not quantitatively very 
accurate in the weak-coupling regime. 



C. Double occupancy 

Finally, we demonstrate that two-particle correlators 
in the charge sector can also be reliably studied with the 
DSR method, taking the fraction of doubly occupied sites 
dfi = {n^ni) as an example. 

Using the constraint (0), we calculate the connected 



charge susceptibility in the following manner: 

Xc(r) = /e(^)-0E(^'(°)-^)} 

= (l(t)L(O)) = ± (d T 0(r)d T e(O)) (41) 
= 772 [Gx(T)d 2 T G x (T) + US(t) - (d T G x (T) f] 

The double occupancy is obtained by taking the equal- 
time value: 

Xc(r = 0) = 2{n t n l )=2d n (42) 

In Fig. jnj we plot the double occupancy obtained in 
that manner, in comparison to the ED and IPT results. 
Within ED, dfj, is calculated directly from the charge 
correlator. Within IPT, Xc( T ) is not approximated very 
reliably Q , but the double occupancy can be accurately 
calculated by taking a derivative of the internal energy 
with respect to U. 

Fig. [ll] demonstrates that the DSR method is very ac- 
curate in the Mott transition region and in the insulator. 
In particular, the hysteretic behaviour is well reproduced. 
In the weak coupling regime however, the approximation 
deteriorates. This issue actually depends on the quan- 
tity: the physical Green's function has the correct limit 
U — > 0, as emphasized above, but this is not true for the 
charge susceptibility (hence for dfjj. The mathemati- 
cal reason is that the constraint (Q) is crucial for writing 
( pd"| ) , but this constraint is only treated on average within 
our method. This shows the inherent limitations of slave- 
boson techniques for evaluating two-particle properties. 
The frequency dependence of two-particle correlators will 
be dealt with in a forthcoming publication. 
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FIG. 11: Double occupancy at f3 = 80 as a function of U. 
The ED result is indicated with squares. DSR overshoots the 
exact result on the left, whereas IPT is slightly displaced to 
the right. 
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FIG. 12: Instability lines U c i(T),U C 2{T), and coexistence re- 
gion for one, two, three and four orbitals (N — 2, 4, 6, 8) as 
obtained by DMFT (DSR solver) 



D. Multi-orbital effects 

In this section, we apply the DSR method to study the 
dependence of the Mott transition on orbital degeneracy 
N. We emphasize that there are at this stage very few 
numerical methods which can reliably handle the multi- 
orbital case, specially when N becomes large. Multi- 
orbital extensions of IPT have been studied || , but the 
results are much less satisfactory than in the one-orbital 
case. The ED method is severely limited by the expo- 
nential growth of the size of the Hilbert space. The Mott 
transition has been studied in the multi-orbital case us- 
ing QMC, with a recent study [||| going up to N = 8 
(4 orbitals with spin). Furthermore, we have recently 
obtained 1 15 some analytical results on the values of the 
critical couplings in the limit of large- N. Those, together 
with the QMC results, can be used as a benchmark of the 
DSR approximation presented here. As explained above, 
a value of Nn=2 — 3 was found to describe best the sin- 
gle orbital case (N = 2). Hence, upon increasing N, we 
choose the parameter Af such that Nn/Nn=2 = N/2. 

In Fig. [l2], we display the coexistence region in the 
(U, T) plane for increasing values of N, as obtained with 
DSR. The values of the critical interactions grow with 
N. A fit of the transition lines U C \(T) (where the insula- 
tor disappears) and U C 2 (T) (where the metal disappears) 
yields: 



We display on Fig. [l3] the DSR result for the spectral 
function for N = 2 and N = 4, at a fixed value of U and 
T. Two main effects should be noted. First, correlations 
effects in the metal become weaker as N is increased (for 
a fixed U), as clear e.g from the increase of the quasi- 
particle bandwith. This is due to increasing orbital fluc- 
tuations. Second, the Hubbard bands also shift towards 
larger energies, an effect which can be understood in the 
way atomic states broadens in the insulator J20| . 

To conclude this section, we have found that the DSR 
yields quite satisfactory results when used in the DMFT 
context for multi-orbital models. In the future, we plan 
to use this method in the context of realistic electronic 
structure calculations combined with DMFT, in situa- 
tions where "numerically exact" solvers (e.g QMC) be- 
come prohibitively heavy. There are however some limi- 
tations to the use of DSR (at least in the present version 
of the approach), which are encountered when the occu- 
pancy is not close to N/2. We examine this issue in the 
next section. 



E. Effects of doping 



U cl (N,T) = A 1 (T)V N 
U c2 (N,T) = A 2 {T)N 



(43) 
(44) 



These results are in good accordance with both the QMC 
data [^H , and the exact results established in jl5) . When 
increasing the number of orbitals, the coexistence region 
widens and the critical temperature associated with the 
endpoint of the Mott transition line also increases. 



Up to now, we studied the half-filled problem (i.e. con- 
taining exactly N/2 electrons), with exact particle-hole 
symmetry. In that case, eo = (d + U/2 = and the La- 
grange multiplier h enforcing the constraint ( |3l| ) could be 
set to h = 0. In realistic cases, particle-hole symmetry 
will be broken (so that eo 7^ 0), and we need to consider 
fillings different from N/2. Hence, the Lagrange multi- 
plier h must be determined in order to fulfill (J31[), which 
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FIG. 13: Spectral function at U = 2, (3 = 60 and for one (full 
line) and two orbitals (broken line) 



we rewrite more explicitly as: 
1 2h 



MV 



(45) 



-iu n ( e w " 0+ + e iv ^~ 



— - V 

+ AfU vl/U + \-2ihv n /U -Y,x{iv n ) 

In this expression, n f = ^ £ = s Ea(44) 

is the average occupancy per orbital flavor. (Note that 
the number of auxiliary fermions and physical fermions 
coincide, as clear from (j34|)). n/ is related to the d-level 
position (or chemical potential, in the DMFT context) 
by: 



away from half-filling. A critical value of the chemical 
potential (i.e. of eo) is required to enter the metallic 
state. The spectral function displays the three expected 
features: a lower and upper Hubbard bands, as well as 
a quasiparticle peak (which in this case where the dop- 
ing is rather large, is located almost at the top of the 
lower Hubbard band). However, it is immediately appar- 
ent from this figure that the DSR method in the doped 
case overestimates the spectral weight of the upper Hub- 
bard band. This can be confirmed by a comparison to 
other solvers (e.g ED). Note that the energy scale be- 
low which a causality violation appears within U-NCA 
becomes rapidly large as the system is doped, while no 
such violation occurs within DSR. 

The DSR method encouters severe limitations however 
as the total occupancy becomes very different from N/2. 
This is best understood by studying the dependence of 
the occupancy upon eo, in the atomic limit. As explained 
in Appendix. |A|, the use of sigma-model variables X (in 
the large- M limit) results in a poor description of the rif 
vs. eo dependence ("Coulomb staircase"). As a result, 
it is not possible to describe the Mott transitions occur- 
ring in the multi-orbital model at integer fillings differ- 
ent from N/2 using the DSR approximation. It should 
be emphasized however that this pathology is only due 
to the approximation of the 0(2) quantum rotor e by 
a 0(M) sigma-model field. A perfect description of the 
Coulomb staircase is found for all fillings when treating 
the constraint (^J) on average while keeping a true quan- 
tum rotor [^7]. Hence, it seems feasible to overcome this 
problem and extend the practical use of the (dynamical) 
slave rotor approach to all fillings. We intend to address 
this issue in a future work. 



e + h - S/(jw n ) 



(46) 



VI. CONCLUSION 




FIG. 14: Doping the insulator (U = 3 and (3 = 60) with 
eo = 0, 1 (corresponding to nf — 0.5, 0.4 respectively) 

We display in Fig. |l4| the spectral function obtained 
with the DSR solver when doping the Mott insulator 



We conclude by summarizing the strong points as well 
as the limitations of the new quantum impurity solver 
introduced in this paper, as well as possible extensions 
and applications. 

On the positive side, the DSR method provides an 
interpolating scheme between the weak coupling and 
atomic limits (at half- filling) . It is also free of some of 
the pathologies encountered in the simplest finite-?/ ex- 
tensions of NCA (negative lifetimes at low temperature). 
When applied in the context of DMFT, it is able to re- 
produce many of the qualitatively important features as- 
sociated with the Mott transition, such as coexisting in- 
sulating and metallic solutions and the existence of two 
energy scales in the DMFT description of a correlated 
metal (the quasi-particle coherence bandwith and the 
"preformed gap"). Hence the DSR solver is quite useful 
in the DMFT context, at a low computational cost, and 
might be applicable to electronic structure calculations 
for systems close to half-filling when the orbital degen- 
eracy becomes large. To incorporate more realistic mod- 
elling, one can introduce different energy levels for each 
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correlated orbital, while the extension to non-symmetric 
Coulomb interactions (such as the Hund's coupling) may 
require some additional work. 

The DSR method does not reproduce Fermi-liquid be- 
haviour at low energy however, which makes it inad- 
equate to address physical properties in the very low- 
energy regime (as is also the case with NCA). The main 
limitation however is encountered when departing from 
half-filling (i.c from N/2 electrons in an TV-fold degener- 
ate orbital). While the DSR approximation can be used 
at small dopings, it fails to reproduce the correct atomic 
limit when the occupancy differs significantly from N/2 
(and in particular cannot deal with the Mott transition 
at other integer fillings in the multi-orbital case). We 
would like to emphasize however that this results from 
extending the slave rotor variable to a field with a large 
number of components. It is possible to improve this 
feature of the DSR method by dealing directly with an 
0(2) phase variable, which does reproduce accurately the 
atomic limit even when the constraint is treated at the 
mean-field level. We intend to address this issue in a 
future work. Another possible direction is to examine 
systematic corrections beyond the saddle-point approxi- 
mation in the large N, M expansion. 

Finally, we would like to outline some other possible 
applications of the slave rotor representation introduced 
in this paper (Sec. II). This representation is both phys- 
ically natural and economical. In systems with strong 
Coulomb interactions, the phase variable dual to the lo- 
cal charge is an important collective field. Promoting 
this single field to the status of a slave particle avoids 
the redundancies of usual slave-boson representations. 
In forthcoming publications, we intend to use this rep- 
resentation for: i) constructing impurity solvers in the 
context of extended DMFT Q , in which the frequency- 
dependent charge correlation function must be calculated 
p8j ii) constructing mean-field theories of lattice models 
of correlated electrons (e.g the Hubbard model) and 
iii) dealing with quantum effects on the Coulomb block- 
ade in mesoscopic systems. 



perature (at finite temperature, deviations from the exact 
result are of order exp (—/3U) and therefore negligeable 
for pratical purposes). To do this, we first extract the val- 
ues of the mean-field parameters A and h from the saddle- 
point equations at zero temperature and A(r) = 0. 



1 = 



1 



2vr 



+ A + 



2ihv 



(Al) 



2h 



4/? 



AfU 2 / 2tt fid 



Performing the integrals shows that A = (U 2 — 4ft, 2 )/(4£7) 
and 6(h - e ) = 1/2, so that h = e . If |e | > U/2, the 
equations lead actually to a solution with an empty or 
full valence that we show on Fig. [l5| . 

We can now compute the physical Green's function 
Gd{T) = G f{r)Gx{—T) from the pseudo-propagators 
Gf(iuj n ) = l/(ito n ) and: 



G x {iv n ) = 



1 



vl/XJ + (U/4 - e 2 Q /U) - 2ie Q v n /U 
-1 1 



(A2) 



iv n + £q - U/2 iv n + e + U/2 



Performing the convolution in imaginary frequency 
and taking the limit T = leads to: 



G d (iuj n ) = -^y^^G x {iv n )Gf{ 



iuj n + iu n . 



1/2 



1/2 



iuj n — eo + U/2 iuj n — 6q — U/2 



(A3) 



(A4) 



Because eo = —fi + U/2 this is the correct atomic limit of 
the single-band model (at half- filling) . The result for the 
empty or full orbital is however not accurate, as shown in 
figure [H| This discrepancy with the correct result (even 
for one orbital) finds its root in the large M treatment 
of the slave rotor X. 
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APPENDIX A: THE ATOMIC LIMIT 

In this appendix we prove the claim that the atomic 
limit of the model is exact at half-filling and at zero tem- 



APPENDIX B: NUMERICAL SOLUTION OF 
THE REAL TIME EQUATIONS 



Here we show how the saddle-point equations (£8[]31j) 
can be analytically continued along the real axis. We 
start with E^(t) = — A/"A(r)G/(— r), which can be first 
Fourier transformed into: 



Jo 



(Bl) 



dei 



de 



= - N / ^G';( ei ) / ^A"(e 2 ) 

7 7T J 7T W n + €1 - £ 2 

where we used the spectral representation 

duj G"(e) 



n F (ei) - n F {e 2 ) 



G{z) 



it z — e 



(B2) 
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FIG. 15: Impurity occupancy in function of the d-level posi- 
tion in the rotor description (full curve) and the exact result 
(dot curve), in the two-orbital case. 



where n/, the average number of physical fermions, is: 

n f = G/(r = 0-) = -J ^G'}(e)n F (e) (B9) 

We note here that solving these real-time integral equa- 
tions can be quite difficult deep in the Kondo regime of 
the Anderson model, or very close to the Mott transition 
for the full DMFT equations. The reason is that Gf(uj) 
and G x (tu) develop low energy singularities (this is an- 
alytically shown in the next appendix), that make the 
numerical resolution very unprecise if one uses FFT. In 
that case, it is necessary to introduce a logarithmic mesh 
of frequency (loosing the benefit of the FFT speed, but 
increasing the accuracy), or to perform a Pade extrapo- 
lation of the imaginary time solution. 



APPENDIX C: FRIEDEL'S SUM RULE 



for each Green's function. The notations here are quite 
standard: G"(e) = TmG(e + i0 + ), n F (e) is the Fermi 
factor, and ng(e) denotes the Bose factor. 

It is then immediate to continue iv n — > v+iO + in equa- 



tion (Bl), and using again the spectral decomposition 
(B2), we derive an equation between retarded quantities: 



G'}(e)n F (e)A(e + v) (B3) 



A" / —A"(e)n F (e)G f (e-v) 

7T 



A calculation along the same lines for the fermionic self- 
energy £/(t) = A(t)G x (t) leads to: 

E/M = -/ ^G' x (e)n B (e)A^-e) (B4) 



Ac 



A"(e)n F (e)G x (cj-e) 



The n ume rical i mpl ementation is then straightforward 
because (B3) and (B4) can each be expressed as the con- 
volution product of two quantities, so that they can be 
calculated rapidly using FFT. The algorithm is looped 
back using the Dyson equations (for real frequency): 



G f 1 (w) = uj — e + h — 



G x \u) 



U 



A 



2hv 

IT 



(B5) 
ExO) (B6) 



At each iteration, A and h are determined using a bi- 
section on the equations ( |30f]3"l| ) , which can be properly 
expressed in terms of retarded Green's functions: 

1 = J ^ G ' x {e)n B {e) ( B7) 

J ^G' x (e)en B (e) (B8) 



2h 

JfU~ATU 



We present here for completeness the derivation of the 
Slave Rotor FriedePs sum rule, equation (|37|), at half- 
filling. The idea, motivated by the numerical analysis as 
well as theoretical arguments [[l2[ |3^| , is that the pseudo- 
particles develop low frequency singularities at zero tem- 
perature: 



G» = Af\w\~ a t 

G' x {uo) = AjtM-^SignCw) 

Using the spectral representation 

G(r) = [ +C ° —e^ T G"{uj) 



(CI) 
(C2) 



(C3) 



we deduce the long time behavior of the Green's functions 
(we denote by T(z) the gamma function): 



Gf(r) 
Gx(r) 
We have similarly 



AjlXl - a/ ) Sign(r) 
7r \r\ 1 - a f 
A x T(l-a x ) 1 



A(r) 



A"(0) 



(C4) 
(C5) 

(C6) 



if one assumes a regular bath density of states at zero 
frequency. The previous expressions allow to extract the 
long time behavior of the pseudo-self-energies (using the 
saddle-point equations (p8f]29|)): 



Ex(r) 
E/(r) 



NA f r(l - of) A"(0) 
it 2 \r\ 2 ' a f 
A X T{1 - a x ) A"(0) Sign(r) 



(C7) 
(C8) 
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The n ext step is to use ( |C3 ) the other way around to get 
from (C7-C8|) the w-depcndance of the self-energies: 



,„,, _ ^^/A"(0) |w|1 _ a/Sigii(a;) (C9) 



7r 1 — a.f 



S/M = - 



UxA"(0). 



7r 1 — ax 



(CIO) 



It is necessary at this point to calculate the real part of 
both self-energies. This can be done using the Kramers- 
Kronig relation, but analyticity provides a simpler route. 
Indeed, by noticing that is an analytic function of 
z and must be uni-valuated above the real axis , we find: 

» M ^/A"(o) ^ 

* l ' ~ 7T l- a/ sin[(a/ - 1)tt/2] 1 ' 

= 1 A X A (0) e u|i-*x 

7r 1 — ax sm[ax7r/2J 

The same argument shows from (|C1-|C2) that: 

i(a / +l)7r/2 

C ' I I — nt 



Gf(z) 



A, 



sm[(a f + l)7r/2] 1 



Gx(z) = Ax 



ixt/2 



sin [ax tt/ 2 



(Cll) 
(C12) 



We can therefore collect the previous expressions, using 
Dyson's formula for complex argument: 



G x \z) = -f-+A-SxW 



U 



(C13) 
(C14) 



and this enables us to extract the leading exponents, as 
well as the product of the amplitudes: 



1 



ax 
A f A x 



TV + 1 

N 
Af+1 



TV + 1 A"(0) 



sin — 



7T TV 



2N+1 



(C15) 
(C16) 
(C17) 



We finish by computing the long time behavior of the 
physical Green's function GA t) = G/(r)Gx(— t) to- 
gether with equations (C4-C5|): 



G *> = WW a " {iTrr^js^rr < C18) 



This proves (f57|). In principle, next leading order cor- 
rections can be computed by the same line of arguments, 
although this is much more involved p2[ . Non Fermi Liq- 
uid correlations in the physical Green's function would 
appear in this computation. 
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